home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / vol_200 / 212_01 / perc3dr.c < prev    next >
Encoding:
C/C++ Source or Header  |  1980-01-01  |  13.5 KB  |  604 lines
  1. /* PERC3DR.C    VERS:- 01.00  DATE:- 10/06/86  TIME:- 05:23:48 PM */
  2. /*
  3. %CC1 $1.C -X -O -E 7000
  4. %CLINK $1 DIO -F FLOAT -S
  5. %DELETE $1.CRL 
  6. */
  7. /* 
  8. Description:
  9.  
  10. Simulation of percolation on a three-dimensional simple cubic lattice.
  11. Option to wrap in one dimension, ie to simulate cylindrical shells.
  12. The algorithm is fairly efficient (20 seconds/1000 site 2-d sheet).
  13.  
  14. A percolative process is characterized by a percolation transition, a phase
  15. change at which long-range connectivity (e.g., a gel phase; conduction 
  16. across the sample; etc) suddenly appears.  Connectivity is defined by the
  17. extent of a cluster of filled sites adjacent (ie, connected) in the lattice.
  18. Long-range connectivity is equivalent to a cluster of infinite extent.
  19.  
  20. The percolation model applies to stochastic processes taking place 
  21. in a spatially disordered system.  
  22. Examples: 
  23.     gelation; 
  24.     flow in a porous medium; 
  25.     dilute magnets; 
  26.     conductivity of conductor-insulator composite materials.
  27.  
  28. For a discussion of percolation, see R. Zallen, "The Physics of Amorphous
  29. Solids", John Wiley, New York, 1983.
  30.  
  31. This program finds the fractional occupancy of the model lattice at the 
  32. percolation transition.  Lattice sites are filled randomly.  When each filled
  33. site is added to the lattice, connectivity is tested along one dimension
  34. of the lattice (ie, between one pair of arbitrarily chosen sides of the 
  35. lattice).  The point of the percolation transition is the point at which
  36. connectivity between sides first appears.
  37.  
  38. Because of the stochastic character of percolation, a simulation consists
  39. of the averaging of NRPT determinations of the percolation transition.
  40.  
  41. The program is set for looping over the simulation, with change in
  42. control parameters at each loop.
  43.  
  44. Requires sqrt() function.
  45.  
  46. By J.A. Rupley, Tucson, Arizona
  47. Coded for BDS C compiler, version 1.50a
  48. */
  49.  
  50. #include "jario.h"
  51. #include "bdscio.h"
  52. #include "dio.h"
  53.  
  54. /*CHANGE TO FALSE IF SYSTEM NOT CONFIGURED FOR HAYES MODEM AS PUNCH-READER*/
  55. #define HAYES_MODEM FALSE
  56. /*CHANGE TO TRUE FOR DEBUG OUTPUT*/
  57. #define DEBUG FALSE
  58.  
  59. #define NSITES 6400
  60. #define NROW  10
  61. #define NCOL  10
  62. #define NLEVEL 10
  63. #define NRPT  20
  64. #define NO_PRINT 2
  65. #define NCLUSTERS (NSITES / 8)
  66. #define WRAP_COL FALSE
  67.  
  68. int mast_array[NSITES];
  69. int top[NCLUSTERS], bottom[NCLUSTERS];
  70. int cnt_crit[100];
  71.  
  72. char str_buf[80];
  73. char title[160];
  74.  
  75. int xnrpt, xnrow, xncol, xnlevel;
  76. int nrpt, nrpt_1;
  77. int nsites;
  78. int no_print;
  79. int wrap_col;
  80.  
  81. int crit_sum;
  82. char fcrit_avg[5], fcrit_var[5], fcrit_1var[5];
  83. char fract_avg[5], fract_var[5], fract_1var[5];
  84. char fnrpt[5], fnrpt_1[5], ftemp[5], fsites[5], fsqd[5];
  85.  
  86. int i, j, k, ii, jj, kk, ll;
  87. int i_new, j_new, k_new;
  88. int i_old, j_old, k_old;
  89. int new_index, old_index;
  90. int new, old;
  91. int cnt_cluster, ncluster;
  92. char pspecial;
  93. int s1, s2, s3;
  94. char clock_in[12];
  95.  
  96. /* page eject*/
  97.  
  98.  
  99. main(argc, argv)
  100. int argc;
  101. char **argv;
  102. {
  103.     dioinit(&argc, argv);
  104.     /*PSPECIAL - LOOP OVER VARIED PARMS*/
  105.     /*SETUP*/
  106.     for (pspecial = 0; pspecial < 100; pspecial++)
  107.     {
  108.         switch (pspecial)
  109.         {
  110.         case 0 :        /*case 0 used -ONLY- for setup*/
  111.             setup();
  112.             /*insert title(up to 2 lines)*/
  113.             strcpy(title,
  114.             "SIMULATION OF PERCOLATION ON A THREE-DIMENSIONAL SIMPLE CUBIC LATTICE"
  115.                 );
  116.             strcat(title,
  117.             ""
  118.                 );
  119.             printf("\n%s\n", title);
  120.             read_variables();
  121.             set_more();
  122.             continue;
  123.         case 1 :        /*case 1 to 99 = specials*/
  124.             break;
  125.         default :
  126.             dioflush();
  127.             exit();
  128.         }
  129.  
  130.         header_disp();
  131.  
  132.         /*MAIN LOOP*/
  133.         /*repeat determination of critical count */
  134.         /*to obtain the average*/
  135.         for (nrpt = 0; nrpt < xnrpt; nrpt++)
  136.         {
  137.             /*initialize mast_array */
  138.             /*with optional blocking of sites*/
  139.             init_array();
  140.  
  141.             /*SUB LOOP*/
  142.             /*start of cycle of filling sites,*/
  143.             /*until have connection top-to-bottom*/
  144.             /*the critical count is the number filled*/
  145.             /*at the point of first connection*/
  146.             while (TRUE)
  147.             {
  148.                 /*find unfilled site, by random search*/
  149.                 i_new = (nrand(1) % xnrow);
  150.                 j_new = (nrand(1) % xncol);
  151.                 k_new = (nrand(1) % xnlevel);
  152.                 new_index = i_new * xncol * xnlevel
  153.                     + j_new * xnlevel + k_new;
  154.                 if (mast_array[new_index])
  155.                     continue;
  156.                 /*add new site to the master array*/
  157.                 /*exit subloop if connection made*/
  158.                 if (add_site())
  159.                     break;
  160. #if !DEBUG
  161.                 if (exit_test())
  162.                     array_disp();
  163. #endif
  164.             }
  165.  
  166.             /*display cnt_crit values and their*/
  167.             /*average and variance*/
  168.             if (nrpt < no_print)
  169.             {
  170.                 printf("\n\n");
  171.                 array_disp();
  172.                 crit_disp();
  173.             }
  174.             else
  175.                 typef("working on cycle %d\n", (nrpt + 1));
  176. #if DEBUG
  177.             exit_test();
  178. #endif
  179.         }
  180.         /*ALL DONE--TERMINATION OF MAIN LOOP*/
  181.         prog_exit();
  182.     }
  183.     /*END OF PSPECIAL LOOP*/
  184. }
  185. /* END OF MAIN                */
  186. /* page eject*/
  187.  
  188.  
  189. void read_variables()
  190. {
  191.     printf("\nDo you want to change values of program variables? (y/cr=>default): ");
  192.     getline(str_buf, 2);
  193.     if ((str_buf[0] != 0x79) && (str_buf[0] != 0x59))
  194.     {
  195.         printf("\n");
  196.         return;
  197.     }
  198.     printf("\nNew title? \n*:");
  199.     getline(str_buf, 79);
  200.     if (strlen(str_buf))
  201.     {
  202.         strcpy(title, str_buf);
  203.         printf("     Another title line?\n*:");
  204.         getline(str_buf, 79);
  205.         if (strlen(str_buf))
  206.             strcat(title, "\n");
  207.         strcat(title, str_buf);
  208.     }
  209.     printf("\n");
  210.     readval("How many array rows? (value/cr=>no change = %d): ", &xnrow);
  211.     printf("\n");
  212.     readval("How many array cols? (value/cr=>no change = %d): ", &xncol);
  213.     printf("\n");
  214.     readval("How many array levels? (value/cr=>no change = %d): ", &xnlevel);
  215.     printf("\n");
  216.     readval("Wrap left column to right column (1 = yes)? (value/cr=>no change = %d): ", &wrap_col);
  217.     printf("\n");
  218.     readval("How many repeat determinations? (value/cr=>no change = %d): ", &xnrpt);
  219.     printf("\n");
  220.     readval("How many cycles do you want printed? (value/cr=>no change = %d): ", &no_print);
  221.     printf("\n\n");
  222. }
  223. /*END OF READ_VARIABLES            */
  224.  
  225.  
  226.  
  227. void readval(string, val)
  228. char *string;
  229. int *val;
  230. {
  231.     while (TRUE)
  232.     {
  233.         printf(string, *val);
  234.         if (strlen(gets(str_buf)) == 0)
  235.             return;
  236.         *val = atoi(str_buf);
  237.     }
  238. }
  239. /*END OF READVAL            */
  240.  
  241.  
  242.  
  243. void setup()
  244. {
  245.     /* set default values of xnrow, etc*/
  246.     xnrow = NROW;
  247.     xncol = NCOL;
  248.     xnlevel = NLEVEL;
  249.     xnrpt = NRPT;
  250.     no_print = NO_PRINT;
  251.     wrap_col = WRAP_COL;
  252.  
  253.     /*establish kernel for random number generator*/
  254.     read_clock();
  255.     nrand(-1, s1, s2, s3);
  256. }
  257. /*END OF SETUP                */
  258.  
  259.  
  260.  
  261. void set_more()
  262. {
  263.     if ((nsites = xnrow * xncol * xnlevel) > NSITES)
  264.     {
  265.         printf("\n\nArray too large.  Error exit.");
  266.         dioflush();
  267.         exit();
  268.     }
  269. }
  270. /*END OF SET_MORE                */
  271.  
  272.  
  273. #if HAYES_MODEM
  274. void read_clock()
  275. {
  276.     /* Clock_in gives hhmmssx    */
  277.     hayes_chrono("RT\r", clock_in);
  278.     s1 = atoi(&clock_in[2]);
  279.     printf("\nTime, date and day =  %s ", clock_in);
  280.  
  281.     /* Clock_in gives yymmdd    */
  282.     hayes_chrono("RD\r", clock_in);
  283.     s2 = atoi(&clock_in[2]);
  284.     printf("%s ", clock_in);
  285.  
  286.     /* Clock_in gives  d = acsii #    */
  287.     hayes_chrono("RW\r", clock_in);
  288.     s3 = atoi(clock_in);
  289.     printf("%s \n", clock_in);
  290. }
  291. /*END OF READ_CLOCK            */
  292.  
  293.  
  294.  
  295. #define PUNCH 6
  296. #define READER 7
  297. #define CR 0x0d
  298.  
  299. void hayes_chrono(command, buffer)
  300. /* Send command to clock and return with string*/
  301. char *command, *buffer;
  302. {
  303.     strcpy(str_buf, "AT");
  304.     strcat(str_buf, command);
  305.     for (i = 0; str_buf[i] != NULL; bios(PUNCH, str_buf[i++]))
  306.         ;
  307.     for (i = 0; i < 12; i++)
  308.         if ((str_buf[i] = bios(READER, 0)) == CR)
  309.             break;
  310.     str_buf[i] = '\0';
  311.     strcpy(buffer, str_buf);
  312. }
  313. /* END OF HAYES_CHRONO    */
  314.  
  315.  
  316.  
  317. #else
  318. void read_clock()
  319. {
  320.     printf("\n\nPlease hit three keys, not too quickly, to seed random number generator\n\n");
  321.     s1 = s2 = s3 = 0;
  322.     while (!bios(2, 0))
  323.         s1++;
  324.     bios(3, 0);
  325.     while (!bios(2, 0))
  326.         s2++;
  327.     bios(3, 0);
  328.     while (!bios(2, 0))
  329.         s3++;
  330.     while (bios(2, 0))
  331.         bios(3, 0);
  332. }
  333. /*END OF "READ_CLOCK" FOR SYSTEM WITH NO CLOCK */
  334. #endif
  335.  
  336.  
  337.  
  338. int exit_test()
  339. {
  340.     int t;
  341.  
  342. #if !DEBUG
  343.     if (bios(2, 0))
  344. #endif
  345.     {
  346.         while (bios(2, 0))
  347.             bios(3, 0);
  348.         printf("\ntype x or X to exit or any other character to continue: ");
  349.         if (((t = bios(3, 0)) == 0x58) || (t == 0x78))
  350.         {
  351.             prog_exit();
  352.             dioflush();
  353.             exit();
  354.         }
  355.         return TRUE;
  356.     }
  357.     return FALSE;
  358. }
  359. /*END OF EXIT_TEST            */
  360.  
  361.  
  362.  
  363. void prog_exit()
  364. {
  365.     header_disp();
  366.     nrpt = nrpt - 1;
  367.     crit_disp();
  368.     printf("\n\n\n\f");
  369. }
  370. /*END OF PROG_EXIT            */
  371.  
  372.  
  373.  
  374. void header_disp()
  375. {
  376.     printf("\n\n%s\n", title);
  377. #if HAYES_MODEM
  378.     read_clock();
  379. #endif
  380.     printf("\nSimulation of percolation on a three-dimensional simple cubic lattice.");
  381.     printf("\nThe array size is %d rows by %d columns by %d levels, %d sites total.",
  382.         xnrow, xncol, xnlevel, nsites);
  383.     printf("\nThe test is for connection between top and bottom levels.");
  384.     if (wrap_col)
  385.     {
  386.         printf("\nWrap left to right column.");
  387.         printf("\nCylindrical shell: diameter = columns; thickness = rows; height = levels.\n\n");
  388.     }
  389.     else
  390.         {
  391.         printf("\nWrap off.");
  392.         printf("\nPlate: width = columns; thickness = rows; height = levels.\n\n");
  393.     }
  394. }
  395. /*END OF HEADER_DISP            */
  396.  
  397.  
  398.  
  399. void crit_disp()
  400. /*display critical counts and fractions--averages and variances*/
  401. {
  402.     char fpno1[5], fpno2[5];
  403.     int sign;
  404.     char *sqrt();
  405.  
  406.     crit_calc();
  407.  
  408.     printf("\nList of all counts of filled sites at the critical point\n");
  409.     for (i = 0; i < (nrpt + 1); i++)
  410.         printf("%4d%c", cnt_crit[i], ((i % 10) == 9 || i == nrpt) ? '\n' : ' ');
  411.     printf("\n\n   critical count =%10.2f     sample_sd=%10.2f   avg_sd=%10.2f", fcrit_avg, sqrt(fpno1, &sign, fcrit_var), sqrt(fpno2, &sign, fcrit_1var));
  412.     printf("\n   critical/total =%10.5f     sample_sd=%10.7f   avg_sd=%10.7f", fract_avg, sqrt(fpno1, &sign, fract_var), sqrt(fpno2, &sign, fract_1var));
  413.     printf("\n\n\n");
  414. }
  415. /* END OF CRIT_DISP            */
  416.  
  417.  
  418.  
  419. void crit_calc()
  420. /*calculate averages and variances of critical counts and fractions*/
  421. {
  422.     nrpt_1 = nrpt + 1;
  423.     itof(fnrpt, nrpt);
  424.     itof(fnrpt_1, nrpt_1);
  425.     itof(fsites, nsites);
  426.     fpmult(fsqd, fsites, fsites);
  427.  
  428.     crit_sum = 0;
  429.     for (i = 0; i < nrpt + 1; i++)
  430.         crit_sum = crit_sum + cnt_crit[i];
  431.  
  432.     itof(fcrit_avg, crit_sum);
  433.     fpdiv(fcrit_avg, fcrit_avg, fnrpt_1);
  434.     atof(fcrit_var, "0");
  435.     if (nrpt > 0)
  436.         for (i = 0; i < (nrpt + 1); i++)
  437.         {
  438.             itof(ftemp, cnt_crit[i]);
  439.             fpsub(ftemp, ftemp, fcrit_avg);
  440.             fpmult(ftemp, ftemp, ftemp);
  441.             fpadd(fcrit_var, ftemp, fcrit_var);
  442.         }
  443.     fpdiv(fcrit_var, fcrit_var, fnrpt);
  444.     fpdiv(fcrit_1var, fcrit_var, fnrpt_1);
  445.     fpdiv(fract_avg, fcrit_avg, fsites);
  446.     fpdiv(fract_var, fcrit_var, fsqd);
  447.     fpdiv(fract_1var, fcrit_1var, fsqd);
  448.  
  449. }
  450. /*END OF CRIT_CALC            */
  451.  
  452.  
  453. void init_array()
  454. {
  455.     /*zero arrays*/
  456.     ll = nsites;
  457.     while (ll--)
  458.         mast_array[ll] = 0;
  459.     cnt_cluster = 0;
  460.     ncluster = 0;
  461.     cnt_crit[nrpt] = 0;
  462. }
  463. /*END OF INIT_ARRAY            */
  464.  
  465.  
  466.  
  467. int add_site()
  468. {
  469.     int wrap_flag;
  470.     cnt_crit[nrpt]++;
  471.     new = 0;
  472.     for (i = 0; i < 6; i++)
  473.     {
  474.         i_old = i_new;
  475.         j_old = j_new;
  476.         k_old = k_new;
  477.         wrap_flag = 0;
  478.         switch (i)
  479.         {
  480.         case 0 :
  481.             i_old = max(0, i_new - 1);
  482.             break;
  483.         case 1 :
  484.             i_old = min(xnrow - 1, i_new + 1);
  485.             break;
  486.         case 2 :
  487.             if ((j_old = (j_new - 1)) < 0)
  488.             {
  489.                 wrap_flag = 1;
  490.                 if (wrap_col)
  491.                     j_old = xncol - 1;
  492.                 else
  493.                     j_old = 0;
  494.             }
  495.             break;
  496.         case 3 :
  497.             if ((j_old = (j_new + 1)) >= xncol)
  498.             {
  499.                 wrap_flag = 1;
  500.                 if (wrap_col)
  501.                     j_old = 0;
  502.                 else
  503.                     j_old = xncol - 1;
  504.             }
  505.             break;
  506.         case 4 :
  507.             k_old = max(0, k_new - 1);
  508.             break;
  509.         case 5 :
  510.             k_old = min(xnlevel - 1, k_new + 1);
  511.             break;
  512.         }
  513.         old_index = i_old * xncol * xnlevel
  514.             + j_old * xnlevel + k_old;
  515.         old = mast_array[old_index];
  516. #if DEBUG
  517.         printf("wrap = %d %d old_point= %2d %2d %2d  old_index=%-3d old=%-3d new_index=%-3d new=%-3d\n", wrap_col, wrap_flag, i_old, j_old, k_old, old_index, old, new_index, new);
  518. #endif
  519.         if (!old)
  520.             continue;
  521.         if (!new)
  522.             new = mast_array[new_index] = old;
  523.         else
  524.             if (new != old)
  525.                 combine_cluster();
  526.     }
  527.     if (!new)
  528.     {
  529.         if (ncluster++ == NCLUSTERS)
  530.         {
  531.             printf("\n\nToo many clusters.  Error exit\n");
  532.             prog_exit();
  533.             dioflush();
  534.             exit();
  535.         }
  536.         cnt_cluster++;
  537.         new = mast_array[new_index] = ncluster;
  538.         top[new] = bottom[new] = 0;
  539.     }
  540.     if (k_new == 0)
  541.         top[new] = 1;
  542.     if (k_new == (xnlevel - 1))
  543.         bottom[new] = 1;
  544. #if DEBUG
  545.     array_disp();
  546. #endif
  547.     if (top[new] && bottom[new])
  548.         return 1;
  549.     else
  550.         return 0;
  551. }
  552. /*END OF ADD_SITE            */
  553.  
  554.  
  555.  
  556. void combine_cluster()
  557. {
  558.     cnt_cluster--;
  559.     ll = nsites;
  560.     while (ll--)
  561.         if (mast_array[ll] != old)
  562.             continue;
  563.     else
  564.         mast_array[ll] = new;
  565.     if (top[old])
  566.         top[new] = 1;
  567.     if (bottom[old])
  568.         bottom[new] = 1;
  569. }
  570. /*END OF COMBINE_CLUSTER        */
  571.  
  572.  
  573.  
  574. void array_disp()
  575. {
  576.     int temp_index;
  577.     char temp;
  578.  
  579.     printf("\ncoordinates of last point = %d %d %d    point = %d of %d sites.\n", i_new, j_new, k_new, cnt_crit[nrpt], nsites);
  580.     printf("cluster = %d  of  (%d - %d) clusters.\n", new, ncluster, (ncluster - cnt_cluster));
  581.     printf("topflag = %d, bottomflag = %d        wrap_col = %d (1 => true)\n", top[new], bottom[new], wrap_col);
  582.     if (((3 + xncol) * xnlevel) >= 79)
  583.         return;
  584.     for (i = 0; i < xnrow; i++)
  585.     {
  586.         printf("\n");
  587.         for (k = 0; k < xnlevel; k++)
  588.         {
  589.             printf(" | ");
  590.             for (j = 0; j < xncol; j++)
  591.             {
  592.                 temp_index = 
  593.                     i * xncol * xnlevel + j * xnlevel + k;
  594.                 temp = 0x40 + mast_array[temp_index];
  595.                 printf("%c", (temp == 0x40 ? '-' : temp));
  596.             }
  597.         }
  598.     }
  599.     printf("\n\n");
  600. }
  601. /*END OF ARRAY_DISP            */
  602.  
  603.  
  604.